W niniejszym raporcie przedstawiono działania podjęte w celu określenia głównych przyczyn zmniejszania się długości śledzi oceanicznych wyławianych w Europie. Predykcja została dokonana w oparciu o model przygotowany na zbiorze danych opisanym w sekcji Zbiór danych. W oparciu o przyjęte miary najlepszy okazał się model z algorytmem Random Forest.
W wyniku przeprowadzonej analizy jako czynniki mające wpływ na zmniejszenie się długości śledzi wskazano zmianę temperatury przy powierzchni wody, a także dostępność planktonu Calanus finmarchicus gat. 1 oraz widłonogów gat. 1.
Do przeprowadzenia analizy oraz wygenerowania raporu zostały wykorzystane następujące biblioteki:
library('knitr')
library('dplyr')
library('corrplot')
library('ggplot2')
library('randomForest')
library('caret')
Celem zapewnienia powtarzalności oraz odtwarzalności zaprezentowanych wyników analizy ustawiono stałą wartość ziarna generatora (ang. seed). Fragment kodu odpowiedzialny za tę funkcjonalność został przedstawiony poniżej.
set.seed(123)
Dane wykorzystane do analizy stanowią pomiary śledzi oraz warunków w jakich żyją z ostatnich 60 lat. Dane te zostały zebrane podczas połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.
Zbiór danych został pobrany ze strony http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/sledzie.csv (dostęp 10-11-2016).
W analizowanym zbiorze danych znajdują się wartości puste, nieznane (NA, ang. Not Avaible), które zostały w tym przypadku oznaczone symbolem ?.
herr_data <-read.csv("sledzie.csv", na.strings = "?", comment.char = "")
Znaczenie kolejnych kolumn reprezentujących atrybuty zostały przedstawione w poniższej tabeli.
| atrybut | informacja | jednostka |
|---|---|---|
| length | długość złowionego śledzia | cm |
| cfin1 | dostępność planktonu Calanus finmarchicus gat. 1 | zagęszczenie |
| cfin2 | dostępność planktonu Calanus finmarchicus gat. 2 | zagęszczenie |
| chel1 | dostępność planktonu Calanus helgolandicus gat. 1 | zagęszczenie |
| chel2 | dostępność planktonu Calanus helgolandicus gat. 2 | zagęszczenie |
| lcop1 | dostępność planktonu widłonogów gat. 1 | zagęszczenie |
| lcop2 | dostępność planktonu widłonogów gat. 2 | zagęszczenie |
| fbar | natężenie połowów w regionie | ułamek pozostawionego narybku |
| recr | roczny narybek | liczba śledzi |
| cumf | łączne roczne natężenie połowów w regionie | ułamek pozostawionego narybku |
| totaln | łączna liczba ryb złowionych w ramach połowu | liczba śledzi |
| sst | temperatura przy powierzchni wody | °C |
| sal | poziom zasolenia wody | Knudsen ppt |
| xmonth | miesiąc połowu | numer miesiąca |
| nao | oscylacja północnoatlantycka | mb |
Dane są zorganizowane wierszowo - każdy wiersz odpowiada pojedynczemu połowowi, natomiast kolumny odpowiadają zarejestrowanym zmiennym (atrybutom).
Zbiór danych składa się z 16 atrybutów oraz 52582 rekordów.
kable(summary(herr_data))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 0 | Min. :19.0 | Min. : 0.0000 | Min. : 0.0000 | Min. : 0.000 | Min. : 5.238 | Min. : 0.3074 | Min. : 7.849 | Min. :0.0680 | Min. : 140515 | Min. :0.06833 | Min. : 144137 | Min. :12.77 | Min. :35.40 | Min. : 1.000 | Min. :-4.89000 | |
| 1st Qu.:13145 | 1st Qu.:24.0 | 1st Qu.: 0.0000 | 1st Qu.: 0.2778 | 1st Qu.: 2.469 | 1st Qu.:13.427 | 1st Qu.: 2.5479 | 1st Qu.:17.808 | 1st Qu.:0.2270 | 1st Qu.: 360061 | 1st Qu.:0.14809 | 1st Qu.: 306068 | 1st Qu.:13.60 | 1st Qu.:35.51 | 1st Qu.: 5.000 | 1st Qu.:-1.89000 | |
| Median :26291 | Median :25.5 | Median : 0.1111 | Median : 0.7012 | Median : 5.750 | Median :21.673 | Median : 7.0000 | Median :24.859 | Median :0.3320 | Median : 421391 | Median :0.23191 | Median : 539558 | Median :13.86 | Median :35.51 | Median : 8.000 | Median : 0.20000 | |
| Mean :26291 | Mean :25.3 | Mean : 0.4458 | Mean : 2.0248 | Mean :10.006 | Mean :21.221 | Mean : 12.8108 | Mean :28.419 | Mean :0.3304 | Mean : 520367 | Mean :0.22981 | Mean : 514973 | Mean :13.87 | Mean :35.51 | Mean : 7.258 | Mean :-0.09236 | |
| 3rd Qu.:39436 | 3rd Qu.:26.5 | 3rd Qu.: 0.3333 | 3rd Qu.: 1.7936 | 3rd Qu.:11.500 | 3rd Qu.:27.193 | 3rd Qu.: 21.2315 | 3rd Qu.:37.232 | 3rd Qu.:0.4560 | 3rd Qu.: 724151 | 3rd Qu.:0.29803 | 3rd Qu.: 730351 | 3rd Qu.:14.16 | 3rd Qu.:35.52 | 3rd Qu.: 9.000 | 3rd Qu.: 1.63000 | |
| Max. :52581 | Max. :32.5 | Max. :37.6667 | Max. :19.3958 | Max. :75.000 | Max. :57.706 | Max. :115.5833 | Max. :68.736 | Max. :0.8490 | Max. :1565890 | Max. :0.39801 | Max. :1015595 | Max. :14.73 | Max. :35.61 | Max. :12.000 | Max. : 5.08000 | |
| NA | NA | NA’s :1581 | NA’s :1536 | NA’s :1555 | NA’s :1556 | NA’s :1653 | NA’s :1591 | NA | NA | NA | NA | NA’s :1584 | NA | NA | NA |
Kolejny etap analizy stanowi oczyszczenie danych, którego istotnym elementem jest przetworzenie brakujących danych. Po wnikliwej analizie surowego zbioru danych zauważono, że brakujące wartości można uzupełnić na podstawie grup utworzonych na podstawie wartości atrybutów totaln oraz recr. W ramach tak zdefiniowanych grup wartości posczególnych atrybutów z brakującymi wartościami są takie same. Aby dodatkowo uniezależnić się od wpływu ewentualnych wartości odstających w ramach grupy (ang. outliers) brakujące wartości są estymowane za pomocą mediany.
na_cols <- names(which(colSums(is.na(herr_data))>0))
fill_na <- function(x) ifelse(is.na(x), median(x, na.rm=TRUE), x)
herr_data_nona <- herr_data %>%
group_by(totaln, recr) %>%
mutate_each(funs(fill_na), one_of(na_cols)) %>%
ungroup
Następną fazę przetwarzania stanowiła analiza korelacji pomiędzy zmiennymi. Wartość współczynników korealcji pomiędzy poszczególnymi atrybutmi został zwizualizowany na poniższym rysunku.
herr_mcorr <- cor(herr_data_nona)
corrplot(herr_mcorr, type = "upper", method="color", addCoef.col = "black", tl.col="black", number.digits=2, number.cex=0.75)
Z analizy powyższej macierzy korelacji wynika, że zachodzi duża korelacja między zmiennymi chel1 oraz lcop1 (0.95), chel2 oraz lcop2 (0.89), a także fbar oraz cumf (0.82). Wobec występowania wysokiej korelacji pomiędzy wspomnianymi zmiennymi atrybuty lcop1, chel2 oraz cumf zostały wykluczone, jako redundante, nie wnoszące nowych informacji do analizy.
herr_data_nona <- herr_data_nona[ , -which(names(herr_data_nona) %in% c("lcop1", "lcop2", "cumf", "xmonth"))]
for(i in 2:5) {
print(
ggplot(herr_data_nona, aes(x = herr_data_nona[i])) +
geom_histogram(bins = 30, fill="#56B4E9") +
labs(x = colnames(herr_data_nona)[i], y = "liczba wartości") +
ggtitle(paste("Histogram atrybutu", colnames(herr_data_nona)[i])) +
theme_light()
)
}
Poniższy wykres przedstawia zmianę długości śledzi w czasie (zgodnie z opisem dane były uporządkowane chronologicznie). Niestety w samym zbiorze, jak i jego opisie nie zostały wyspecyfikowane dokładne ramy czasowe pomiarów.
library("plotly")
plot <- ggplot(herr_data_nona, aes(x=X, y=length)) +
geom_point() +
geom_smooth(method="auto", se=TRUE, color="red") +
labs(title="Wykres zmiany długości śliedzi w czasie", x="Czas", y="Długość śledzia")
ggplotly( plot )
Wykres potwierdza postawiony problem badawczy malejącej od pewnego momentu długości śledzia.
W celu stworzenia i oceny regresora konieczne jest podzielenie zbioru danych na dane uczące, walidujące i testowe. Dane uczące i walidacyjne zostały wykorzystane do budowy i optymalizacji parametrów modelu. Natomiast dane testowe posłużyły do wyboru modelu oraz jego końcowej oceny.
in_training <-
createDataPartition(
# atrybut do stratyfikacji
y = herr_data_nona$length,
# procent w zbiorze uczącym
p = .75,
# zwracamy indeksy
list = FALSE)
training <- herr_data_nona[ in_training,]
testing <- herr_data_nona[-in_training,]
W wyniku przeprowadzonych badań model klasyfikacyjny został utworzony zgodnie z algorytmem Random Forest. Do budowy i opytymalizacji modelu wykorzystano powtórzoną ocenę krzyżową z podziałem na 2 części i 4 powtórzeniami.
regr_ctrl <- trainControl(
# powtórzona ocena krzyżowa
method = "repeatedcv",
# liczba podziałów
number = 2,
# liczba powtórzeń
repeats = 4)
rf_regr <- train(length ~ .,
data = training,
# Algorytm Random Forest
method = "rf",
trControl = regr_ctrl,
ntree = 10,
importance = TRUE)
Ostatecznie, w procesie uczenia jako optymalne zostały dobrane następujące wartości parametrów regresora.
## Random Forest
##
## 39438 samples
## 11 predictor
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 4 times)
## Summary of sample sizes: 19718, 19720, 19720, 19718, 19720, 19718, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared
## 2 1.155766 0.5127557
## 6 1.123606 0.5393820
## 11 1.204624 0.4973097
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 6.
Oceny modelu dokonano na podstawie wcześniej przygotowanego zbioru testowego. Jako miary oceny przyjęto R2 oraz RMSE.
predicted <- predict(rf_regr, testing)
observed <- testing[, 'length']
SS_res <- sum( (observed - predicted)^2 )
SS_tot <- sum( (observed - mean(observed))^2 )
R_squared <- 1 - SS_res/SS_tot
RMSE <- sqrt( mean((observed - predicted)^2) / length(predicted) )
length(predicted)
## [1] 13144
| miara | wartość |
|---|---|
| R2 | NA |
| RMSE | 0.01 |
Ocena ważności atrybutów ma za zadanie wskazać jakie atrybuty (czynniki) miały wpływ na to, że rozmiar śledzi zaczął w pewnym momencie maleć.
varImp(rf_regr)
## rf variable importance
##
## Overall
## X 100.000
## sst 51.647
## chel1 33.223
## cfin1 28.144
## chel2 23.050
## fbar 22.568
## totaln 16.474
## nao 14.276
## cfin2 6.206
## recr 5.410
## sal 0.000
Z powyższej tabeli wynika, że najistotniejszy jest atrybut sst reprezentujący temperaturę przy powierzchni wody. Nieco mniej ważny jest atrybut chel1 reprezentujący dostępność planktonu Calanus finmarchicus gat. 1. Warto zwrócić uwagę, że atrybut chel1 był silnie skorelowany (i z tego powodu odrzuconym z analizy, jako redundantny) atrybutem lcop1 (0.95). Zatem również dostępność planktonu w postaci widłonogów gat. 1 ma wpływ na zmniejszenie długości śledzi.
ToDo: Remove X parameter from analysis; correct data formating
Należy zauważyć, że istotność żadnego z atrybutów nie przekracza r, co może oznaczać, że także inne atrybuty nie uwzględnione w zbiorze danych mogły mieć równie istotny wpływ na zmniejszenie się długości śledzi. Przykładowe czynniki mogące mieć wpływ to zanieczyszczenie wody, czy dostępność światła.